Evaluate and detrend OISST anomaly by 30yrs mean

R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/oisst_an.m

First calculate 1982-2011 SST monthly mean. To evaluate SST anomaly:

1. SST.1 = Monthly SST - SST_30yr_mean

2. SST.anomaly = SST.1 - Linear trend by monthly series from 1982-202205

3. Land mask from R/01-2_OISST_land_seaice_mask.Rmd

## The same code from R/01-2_OISST_land_seaice_mask.Rmd(and can check that testing figure)
land_mask <- fread("unzip -p ../data_src/mapfiles/sea_icemask_025d.csv.zip")[, .(lon, lat, x, y, landmask)]
setorder(land_mask, -y, x)
land_mask[,ry:=rowid(x)]
yy <- unique(land_mask[,.(y, ry)]) ## check, NOTE that in OISST.nc file, y is from Northest (89.875) to Southest (-89.875)

landm <- dcast(land_mask, x ~ ry, value.var = "landmask")
xx <- as.numeric(landm[,1]$x)
landm <- as.matrix(landm[,-1])

clim_years = seq(1982,2011) #for climatology
climyrs<- clim_years[length(clim_years)] - clim_years[1] + 1 #30 yrs
monstr <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

Recalculate_30yrs = FALSE

if (Recalculate_30yrs) {
  plan(multisession)
  options(future.globals.maxSize= 1048576000)

  stm <- future_lapply(1:12, function(j) {
    monj <- fifelse(j<10, paste0("0",j), paste0(j))
    jstr <- monstr[j]
    for (i in clim_years) {
      x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
      names(x)[1] <- jstr 
      x[[1]][which(landm==1)] <- NA_real_
      ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
      x[[1]][which(ice[[1]]==1)] <- NA_real_
      if (i == clim_years[1]) {
        stmx <- x
        datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
      } else {
        stmx <- c(stmx, x)
        datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
      }
    }
    names(stmx) <- rep(jstr, length(names(stmx)))
    stmx <- merge(stmx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>% 
        aggregate(by=paste0(climyrs, " years"), FUN=mean, na.rm=TRUE)
    names(stmx)[1] <- jstr
    stmx <- stmx %>% select(jstr) %>% adrop
    return (stmx)
  })
}

if (Recalculate_30yrs) {
  # Note that in an earlier version used "GLDASp4 land mask" but current version use "GLDASp5"
  # The differnt version of land mask cause difference in the definition of land areas, so that sea-mask
  # and then also cause difference in climatology mean result
  save(stm, file="../data_src/stats/sst_1982_2011month_climatology_nc.RData")
} else {
  load("../data_src/stats/sst_1982_2011month_climatology_nc.RData")
}

stm

library(grid)

plotx <- vector("list", length = 12)
pstrx <- '';
for (j in 1:12) {
  plotx[[j]] <- gplotx(stm[[j]], monstr[j], returnx = TRUE, minz = -2.5, maxz = 34.5)
  pstrx <- paste0(pstrx, "plotx[[", j, "]],") 
}

lay1 <- rbind(c(1,1,2,2,3,3),
              c(4,4,5,5,6,6),
              c(7,7,8,8,9,9),
              c(10,10,11,11,12,12))
evplot <- paste0('grid.arrange(', pstrx, ' layout_matrix=lay1)')

gx <- eval(parse(text=evplot))

To double check if the plots of climatology are right, use existed data that download from NOAA (NetCDF OISST, 0.5d, monthly mean)

• Metadata: https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html

• https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc (downloaded 202206 and the data is 1981/12 to 2022/05)

if (!require("ncdf4")) install.packages("ncdf4"); library(ncdf4)

sst_mnfile = "../data_src/stats/sst.mnmean.nc"
if (!file.exists(sst_mnfile)) {
  tryCatch({
    curl::curl_download("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc", destfile = sst_mnfile)
  }, error = function(e) paste0(e))
} 

nx0 <- nc_open(sst_mnfile) 
print(nx0) ## sst[lon,lat,time]   (Chunking: [360,180,464])
latn1<- ncvar_get(nx0, "lat") # 89.5 - -89.5
lngn1<- ncvar_get(nx0, "lon") # 0 - 359.5
time<- ncvar_get(nx0, "time") # 1981-12.01 - 2020-07-01, we now check every March, August, and December
date<- time %>%  as.Date(origin="1800-01-01 00:00:0.0") 

maridx <- seq(4, 464, by = 12) #check date[maridx]
augidx <- seq(9, 464, by = 12)
decidx <- seq(13, 464, by = 12)

get_sstmx <- function (midx, minz=-2.5, maxz=34.5) {
  dtx <- as.data.table(ncvar_get(nx0, "sst")[,,midx]) %>% 
    .[,.(sstm=mean(value, na.rm=TRUE)), by=.(V1,V2)] %>%
    .[,`:=`(longitude=fifelse(lngn1[V1]>180, lngn1[V1]-360, lngn1[V1]), latitude=latn1[V2])] %>%
    .[,.(longitude, latitude, sstm)]
  setnames(dtx, 1:2 , c("x", "y"))
  
  gx <- gplotx(dtx, var="sstm", returnx=TRUE, minz=minz, maxz=maxz, no_coord=TRUE, maxlim=180) 
  gx <- gx +
  #gx <- ggplot() + 
  #  geom_tile(data=dt, aes_string(x="longitude", y="latitude", fill="sstm"), alpha=0.8) + 
  #  scale_fill_viridis(limits=c(minsst, maxsst)) +
     geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3) #+
  #  coord_sf() + 
  #  xlim(c(-180, 180)) + ylim(c(-90, 90))
  
  return(gx)  
}

gx3 <- get_sstmx(maridx)
gx8 <- get_sstmx(augidx)
gx12 <-get_sstmx(decidx)

gt3 <- gplotx(as.data.table(stm[[3]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Mar", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)  
gt8 <- gplotx(as.data.table(stm[[8]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Aug", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)  
gt12 <- gplotx(as.data.table(stm[[12]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Dec", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)  

lay2 <- rbind(c(1,1,2,2),
              c(3,3,4,4),
              c(5,5,6,6))

grid.arrange(gx3, gt3, gx8, gt8, gx12, gt12, layout_matrix=lay2)

i <- 2022
j <- 4
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
  
x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(x)[1] <- "anom" 
#dim(x[[1]]) #1440. 720
#dim(landm)  #1440, 720
#dim(ice)    #1440, 720
x[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
x[[1]][which(ice[[1]]==1)] <- NA_real_

#xt <- matrix()
#length(xt) <-1440 * 720
#dim(xt) <- c(1440, 720)
## if one of x and stm is NA, then the substraction is NA (so numbers of NA values increases)
x[[1]] <- x[[1]] - stm[[j]][[1]]

## Just check if plot 2022-04
gx <- gplotx(x, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
names(y)[1] <- "anom" 
gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)

layt <- rbind(c(1,1,2,2))
grid.arrange(gx, gy, layout_matrix=layt)

yrng <- seq(1982,2022)
endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))

NEED_RETEST_PRED <- FALSE
if (NEED_RETEST_PRED) {
  
j=4 #just testing, arbitrary select a month
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in yrng) {
  if (i==curryr & j>(currmo-1)) break
  z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
  if (i==yrng[1]) {
    stdrx <- z
    datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
  } else {
    stdrx <- c(stdrx, z)
    datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
  }
}
names(stdrx) <- rep(jstr, length(names(stdrx)))
xt <- merge(stdrx)
st_set_dimensions(xt, 3, values = as.POSIXct(datex), names = "time")
seqx <- rbind(c(`377`= as.numeric(xt[[1]][377,41,])), ## arbitrary select a point, check which(!is.na(xt[[1]][,41,]))
              c(`1112`= as.numeric(xt[[1]][1112,41,])),
              c(`test`= sort(rnorm(41)) + rnorm(41)))
colnames(seqx) <- gsub("-", "_", datex)
rownames(seqx) = c("0377_41", "1222_41", "test")
yhat <- stats::predict(lm(seqx[3,] ~ seq_along(datex)))
seqx[1, 28:30] <- NA_real_ #insert NA
seqx[2, 32:34] <- NA_real_ #insert NA
predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
haty <- apply(seqx, MARGIN=1, FUN=function(y){ #predy use predict will auto exclude NA, so that out length is shorter
    if (all(is.na(y))) return(y)
    yh <- as.numeric(predy(y)) 
    y[!is.na(y)] <- yh
    return(c(as.numeric(y)))
}) %>% as.data.frame()  
all.equal(haty$test, as.numeric(yhat)) #TRUE

#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)

ydet <- as.numeric(detrend(seqx[3,], 'linear'))
all.equal(ydet, as.numeric(seqx[3,]) - as.numeric(yhat)) #TRUE

#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)

plot_pred <- function(x, y, yhat=NULL) {
  if (is.null(yhat)) {
    yhat <- stats::predict(lm(y ~ x))
  } else {
    yh <- stats::predict(lm(y ~ x))
    if (!isTRUE(all.equal(as.numeric(yh), yhat))) {
      print(paste0("Not equal yhat with lm: ", paste(yh, collapse=",")))
    }
  }
  plot(x=x, y=y)
  abline(lm(y ~ x))
  points(x=x, y=yhat, col="red", pch=4)
  points(x=x, y=c(y-yhat), col="green", pch=13)
}

plot_pred(seq_along(datex), seqx[3,], haty$`test`)
}

Recalculate_anom <- FALSE
if (Recalculate_anom) {
  
for (i in yrng) {
  mmx <- fifelse(i==curryr, currmo-1L, 12L)
  for (j in 1:mmx) {
    monj <- fifelse(j<10, paste0("0",j), paste0(j))
    jstr <- monstr[j]
    filet <- paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc")
    if (!file.exists(filet)) {
      #z<- read_stars(paste0("../data_src/oisst/monthly_anom/", i, monj, "_anom.nc")) #Old vers use anom, not sst, so no subtract 30yr mean
      z <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
      
      names(z)[1] <- "anom" 
      z[[1]][which(landm==1)] <- NA_real_
      ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
      z[[1]][which(ice[[1]]==1)] <- NA_real_
      z[[1]] <- z[[1]] - stm[[j]][[1]] #Old-version without this substraction while using _anom directly
      write_stars(z, filet)
    }
  }
}
}

## Just check if plot 2022-04
if (Recalculate_anom) {

gx2 <- gplotx(z, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
#y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
#names(y)[1] <- "anom" 
#gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)

layt <- rbind(c(1,1,2,2))
grid.arrange(gx2, gy, layout_matrix=layt)
}

## It's hard/slow to detrend each x,y along almost 30yr (30x12) year trend. Must read 1440x720x30x12 at once...
Recalculate_anom <- FALSE
yrs <- yrng[length(yrng)] - yrng[1] + 1

if (Recalculate_anom) {
for (j in 1:12) {
  monj <- fifelse(j<10, paste0("0",j), paste0(j))
  jstr <- monstr[j]
  for (i in yrng) {
    if (i==curryr & j>(currmo-1)) break
    
    z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))

    if (i==yrng[1]) {
      stdrx <- z
      datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
    } else {
      stdrx <- c(stdrx, z)
      datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
    }
  }
  names(stdrx) <- rep(jstr, length(names(stdrx)))
  print(paste0("Now in Year-month: ", i," - ", monj, " and have length stdrx: ", length(datex)))
  
  predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
  trd <- merge(stdrx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>% 
    aggregate(by=paste0(yrs, " years"), FUN = function(y) {
       if (all(is.na(y))) return(list(as.numeric(y)))
       y[!is.na(y)] <- as.numeric(predy(y)) #Seems if return value has equal length, the result will be simplified by aggregate if using as.matrix
       return(list(as.numeric(y))) #Old-version bug: not predy(y), should be y. See previous testing chunk
    })
  tk <- array(unlist(trd[[1]]), dim = c(length(datex),1440,720)) #bug fix: NOT each predict have the same length result (if input has NA)
  
  print(paste0("Predict ok with dim(tk): ", paste0(dim(tk), collapse=",")))
  for (k in seq_along(yrng)) {
    if (yrng[k]==curryr & j>(currmo-1)) break
    
    z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", yrng[k], monj, "_anom.nc"))
    #notna1 <- which(!is.na(z[1][[1]]))
    #notna2 <- which(!is.na(tk[k,,]))
    #notna <- intersect(notna1, notna2)
    #z[[1]][notna] <- z[[1]][notna]- tk[1,,][notna]
    z[[1]] <- z[[1]]- tk[1,,] #Any of z[[1]] or tk[1,,] NA will cause NA here, it makes sense.
    names(z)[1] <- "anom"
    
    filet <- paste0("../data_src/oisst/monthly_anom_detrend/", yrng[k], monj, "_anom.nc")
    write_stars(z, filet)
    print(paste0("Now write Year-month: ", yrng[k], " - ", monj, " ok"))
  }
}
}  

yk <- read_stars("../data_src/oisst/monthly_anom/202006_anom.nc")
#zk <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"))
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "06", "_anom.nc"))

names(yk)[1] <- "anom" 
#names(zk)[1] <- "anom"
names(z)[1] <- "anom"
print(range(na.omit(as.data.table(yk)$anom))) #-5.911333 10.598333
#print(range(na.omit(as.data.table(zk)$anom))) #-5.39709 16.39870   #old -8.227043 11.727673
print(range(na.omit(as.data.table(z)$anom)))  #-7.724047 10.483911 #old -8.361936  6.181935

gyk<- gplotx(yk,"anom", minz = -6, maxz = 11.0, returnx=TRUE)
gzk<- plot_anom(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"), #zk,
                "anom", #minz = -8, maxz = 11.0, 
                returnx=TRUE)  
gz <- gplotx(z, "anom", minz = -6, maxz = 16.5, returnx=TRUE)  

layt <- rbind(c(1,1),c(2,2),c(3,3))
grid.arrange(gyk, gzk, gz, layout_matrix=layt)

#gyk
#gzk
#gz

if (!require("layer")) {
  if (!require("remotes")) install.packages("remotes"); library(remotes)
  remotes::install_github("marcosci/layer")
  library(layer)
}

Reload_sst_nc <- FALSE
if (Reload_sst_nc) {
  sst1 <- read_stars(paste0("../data_src/oisst/monthly_sst/202002_sst.nc"))
  names(sst1)[1] <- "sst"

  sst202002 <-  as.data.table(sst1) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    #.[,.(longitude, latitude, heatwave)] %>%
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  sst2 <- read_stars(paste0("../data_src/oisst/monthly_sst/201002_sst.nc"))
  names(sst2)[1] <- "sst"

  sst201002 <-  as.data.table(sst2) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  sst3 <- read_stars(paste0("../data_src/oisst/monthly_sst/200002_sst.nc"))
  names(sst3)[1] <- "sst"

  sst200002 <-  as.data.table(sst3) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  sst4 <- read_stars(paste0("../data_src/oisst/monthly_sst/199002_sst.nc"))
  names(sst4)[1] <- "sst"

  sst199002 <-  as.data.table(sst4) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))
  
  save(sst199002, sst200002, sst201002, sst202002, file="../data_src/stats/sst_1990_2020_stack.RData")
} else {
  load("../data_src/stats/sst_1990_2020_stack.RData")
}

tl0 <- layer::tilt_map(sst199002)
tl1 <- layer::tilt_map(sst200002, y_shift=75)
tl2 <- layer::tilt_map(sst201002, y_shift=150)
tl3 <- layer::tilt_map(sst202002, y_shift=225)

map_list <- list(tl0, tl1, tl2, tl3)

#palettem: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
layer::plot_tiltedmaps(map_list, palette = "turbo")